AMF Test27
M. Maltz
20180927
Logic
Salvage topsoil was removed from a donor site and delivered to three recipient sites in 2015. We sampled all sites (Donor and recipients) prior to topsoil delivery and sequenced AM Fungi from all sites in 2015 and 2017 from manipulative plots, set up in a randomized block design at each site with: control treatments (no topsoil added), a dusting of topsoil, and three levels of topsoil thicknesses (2“, 4” and 6" thick layers of this topsoil, originating from the donor site)), delivered to all recipient sites).
knitr::opts_chunk$set(cache=TRUE, tidy=TRUE, error = FALSE, eval = TRUE, message = FALSE, warning = FALSE, rows.print=5, cols.min.print=4, fig.width=6, fig.height=4.5)Questions:
**Q1: Propagule determination: Do AM fungal communities resemble the donor site, regardless of where soil was delivered? Is provenance a driver of AM fungal community composition?
**Q2: Environmental filtering: Do AM fungal communities resemble the recipient sites, and are more dissimilar to the AM fungal communities from the donor site?
**Q3: Propagule pressure: Do AM fungal communities resemble the donor site more in higher level topsoil treatments, than they do in the control or dusted sites?
Outputs:
Analysis of similarity:
Determine the similarity of AMF communities at each recipient site to the donor site, from where topsoils originated from. Generate heatmaps with distances, to visually determine - which sites/samples are more similar to each other, and which are more dissimilar. Generate NMDS / Pcoa, plotted with Year as shapes and Sites as colors, and Description as …? Run permanova Do AM fungal communities at each site resemble the recipient sites, and are they more dissimilar to the donor sites?
Functional group richness:
Determine the richness of AMF functional groups in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment)
Taxonomic diversity:
Determine the OTU taxa richness (alpha diversity) of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment). Determine the beta diversity of AMF communities in the donor site, from where topsoils originated from, and in the recipient sites pre-topsoil delivery (pre-treatment), and post-topsoil delivery (post-treatment).
Topsoil level treatments
Propagule pressure: Do AM fungal communities resemble the donor site more in thicker topsoil layer treatments than they do in the control or dusted sites? Are the sites with Treat groups = S, more similar to the donor site than F, and are treatment groups S and F more similar to the donor site than T? Is this relationship clinal, with S being the thickest and most similar to donor site, followed by F, and then by T. TO FIGURE OUT - how to address this question statistically??multiple regression with Treat on the X and similarity to Donor Site on the Y?? A permanova with multiple comparisons?
knitr::opts_chunk$set(cache = TRUE, tidy = TRUE, error = FALSE, eval = TRUE,
message = FALSE, warning = FALSE, rows.print = 5, cols.min.print = 4, fig.width = 6,
fig.height = 4.5)Load Required Packages
library(tidyverse)
library(dplyr) ## for data wrangling - %>% function
library(reshape2) ##melt and cast data
library(ggplot2) # plotting
library(data.table)
library(stringr)
library(tidyr) # 'separate' function
library(readxl) #read xlsx files into r on mac computer
library(vegan) # dissimilarity matrix, permanova functions
library(magrittr)
library(cowplot)
library(formatR)
date <- format(Sys.time(), "%Y%b%d")
date## [1] "2018Sep27"** the following is adapted from 2017-07-07_Chapter_2_SSU_IntialClean **
Import data files
# read in OTU table
ssu <- fread("CSS_table_sorted.txt")
str(ssu)## Classes 'data.table' and 'data.frame': 386 obs. of 63 variables:
## $ #OTU ID : chr "denovo538" "denovo1813859" "denovo1031897" "denovo34" ...
## $ MM933M : num 5.29 4.33 7.32 0 0 ...
## $ MM936M : num 6.27 11.24 3.4 0 0 ...
## $ MM941M : num 4.4 2.95 0 0 0 ...
## $ MM943M : num 5.31 4.04 0 0 0 ...
## $ MM9362CU: num 3.13 3.13 8.77 0 1.56 ...
## $ MMS700 : num 2.5 6.82 7.82 0 2.5 ...
## $ MMS701 : num 1.87 2.66 8.27 2.66 0 ...
## $ MMS702 : num 0 2.79 8.24 2.79 0 ...
## $ MMS703 : num 0 2.58 7.18 3.45 0 ...
## $ MMS704 : num 2.75 0 7.43 9.75 0 ...
## $ MMS705 : num 3.28 0 6.99 0 4.76 ...
## $ MMS706 : num 1.83 2.61 8.32 0 0 ...
## $ MMS707 : num 3.03 7.49 7.43 12.04 0 ...
## $ MMS708 : num 0 3.88 7.64 2.97 0 ...
## $ MMS709 : num 5.69 8.79 7.95 0 1.56 ...
## $ MMS710 : num 3.98 3.07 8.5 2.23 0 ...
## $ MMS711 : num 4.44 1.99 8.14 1.99 0 ...
## $ MMS712 : num 2.23 3.07 5.91 7.39 4.75 ...
## $ MMS713 : num 2.87 2.05 8.23 2.05 0 ...
## $ MMS714 : num 4.03 4.34 6.54 2.28 8.78 ...
## $ MMS715 : num 3.03 6.25 4.71 9.76 3.03 ...
## $ MMS716 : num 4.22 2.43 8.93 10.72 0 ...
## $ MMS717 : num 3.12 5.56 8.23 2.28 0 ...
## $ MMS890 : num 4.73 2.89 9.24 0 0 ...
## $ MMS891 : num 4.35 2.55 8.66 0 0 ...
## $ MMS900 : num 0 2.95 7.12 2.95 0 ...
## $ MMS901 : num 2.56 6.47 6.11 2.56 0 ...
## $ MMS902 : num 4.22 7.82 6.94 7.64 0 ...
## $ MMS903 : num 4.52 6.31 9.14 2.86 0 ...
## $ MMS904 : num 3.7 2.32 9.18 2.32 0 ...
## $ MMS905 : num 2.95 6.57 7.42 0 0 ...
## $ MMS906 : num 3.5 3.5 9.15 0 0 ...
## $ MMS907 : num 3.54 4.88 7.56 4.69 0 ...
## $ MMS908 : num 5.26 6.77 9 3.22 0 ...
## $ MMS909 : num 3.76 0 6.25 0 0 ...
## $ MMS910 : num 3.4 0 7.02 6.86 0 ...
## $ MMS911 : num 0 0 7.01 4.39 5.73 ...
## $ MMS912 : num 3.01 2.18 8.48 0 0 ...
## $ MMS913 : num 3.94 1.95 9.29 1.95 0 ...
## $ MMS914 : num 3.83 0 7.66 2.44 4.53 ...
## $ MMS915 : num 0 3.72 6.21 2.02 0 ...
## $ MMS916 : num 4.04 2.03 7.68 4.69 0 ...
## $ MMS917 : num 2.58 0 0 0 7.61 ...
## $ MMS918 : num 2.78 8.18 0 0 0 ...
## $ MMS919 : num 3.79 2.4 8.47 0 0 ...
## $ MMS920 : num 4.31 2.26 6.97 3.1 0 ...
## $ MMS921 : num 4.16 2.6 7.69 3.68 1.43 ...
## $ MMS922 : num 3.6 1.92 7.99 2.71 1.92 ...
## $ MMS923 : num 3.55 3.03 8.85 3.55 0 ...
## $ MMS924 : num 3.68 0 8.09 2.79 0 ...
## $ MMS925 : num 3.69 2.31 7.78 3.16 0 ...
## $ MMS928 : num 3.4 3.4 7.87 2.53 0 ...
## $ MMS929 : num 4.19 4.89 9.17 1.95 0 ...
## $ MMS930 : num 4.75 3.81 9.89 0 0 ...
## $ MMS931 : num 0 4.61 8.61 0 0 ...
## $ MMS932 : num 4.63 0 9.28 0 0 ...
## $ MMS938 : num 2.5 2.5 8.11 0 0 ...
## $ MMS939 : num 6.1 3.73 8.08 0 0 ...
## $ MMS942 : num 0 0 8.1 2.97 0 ...
## $ MMS943 : num 0 2.91 7.07 4.92 0 ...
## $ MMS944 : num 5.13 3.87 8.23 2.96 0 ...
## $ taxonomy: chr "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "No blast hit" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" ...
## - attr(*, ".internal.selfref")=<externalptr>## mapping data
metassu <- fread("mapSal18927.txt", header = TRUE)
metassu$Year <- as.character(metassu$Year)
metassu$Rep <- as.factor(metassu$Rep)
metassu$Treat <- as.factor(metassu$Treat)
str(metassu)## Classes 'data.table' and 'data.frame': 61 obs. of 13 variables:
## $ #SampleID : int 34 35 36 41 46 42 43 6 37 38 ...
## $ BarcodeSequence : chr "GAGATCGTTGGGTGAT" "GAGATCGTTACTGTGC" "GAGATCGTAGAGTGTG" "ATGCTTGGTGGTAACC" ...
## $ LinkerPrimerSequence: logi NA NA NA NA NA NA ...
## $ Location : chr "B5" "C5" "D5" "A6" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "DS" "DS" "DS" "DS" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 1 2 3 4 5 1 2 1 2 3 ...
## $ HL : num 1.673 1.519 0.669 0.519 0.679 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "Donor" "Donor" "Donor" "Donor" ...
## $ SampleID2 : chr "MMS938" "MMS939" "MMS890" "MMS891" ...
## - attr(*, ".internal.selfref")=<externalptr>Clean SSU headers
colnames(ssu)## [1] "#OTU ID" "MM933M" "MM936M" "MM941M" "MM943M" "MM9362CU"
## [7] "MMS700" "MMS701" "MMS702" "MMS703" "MMS704" "MMS705"
## [13] "MMS706" "MMS707" "MMS708" "MMS709" "MMS710" "MMS711"
## [19] "MMS712" "MMS713" "MMS714" "MMS715" "MMS716" "MMS717"
## [25] "MMS890" "MMS891" "MMS900" "MMS901" "MMS902" "MMS903"
## [31] "MMS904" "MMS905" "MMS906" "MMS907" "MMS908" "MMS909"
## [37] "MMS910" "MMS911" "MMS912" "MMS913" "MMS914" "MMS915"
## [43] "MMS916" "MMS917" "MMS918" "MMS919" "MMS920" "MMS921"
## [49] "MMS922" "MMS923" "MMS924" "MMS925" "MMS928" "MMS929"
## [55] "MMS930" "MMS931" "MMS932" "MMS938" "MMS939" "MMS942"
## [61] "MMS943" "MMS944" "taxonomy"head(ssu$taxonomy)## [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"
## [3] "No blast hit"
## [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"
## [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"# rename columns
names(ssu)[1] <- "ssuotu" #rename first column
names(ssu)[names(ssu) == "taxonomy"] <- "ssutaxonomy" # rename column that is currently called taxonomy
# split taxonomy column
`?`(str_match)
ssu$ssukingdom <- str_match(ssu$ssutaxonomy, "k__(.*?);")[, 2]
ssu$ssuphylum <- str_match(ssu$ssutaxonomy, "p__(.*?);")[, 2]
ssu$ssuclass <- str_match(ssu$ssutaxonomy, "c__(.*?);")[, 2]
ssu$ssuorder <- str_match(ssu$ssutaxonomy, "o__(.*?);")[, 2]
ssu$ssufamily <- str_match(ssu$ssutaxonomy, "f__(.*?);")[, 2]
ssu$ssugenus <- str_match(ssu$ssutaxonomy, "g__(.*?);")[, 2]
ssu$ssuspecies <- str_match(ssu$ssutaxonomy, "s__(.*?)$")[, 2]
unique(ssu$ssuspecies)## [1] "NES27"
## [2] "MH3_B"
## [3] NA
## [4] "sp"
## [5] "Acau16"
## [6] "MO_G47"
## [7] "laccatum"
## [8] "acnaGlo2"
## [9] "MO_Ar1"
## [10] "leptoticha"
## [11] "Whitfield_type_7"
## [12] "Torrecillas12b_Glo_G5"
## [13] "Sanchez_Castro12b_GLO12"
## [14] "Winther07_D"
## [15] "MO_G18"
## [16] "intraradices"
## [17] "Alguacil12a_Para_1"
## [18] "Douhan9"
## [19] "Alguacil12b_GLO_G3"
## [20] "Glo39"
## [21] "Alguacil10_Glo1"
## [22] "Ligrone07_sp"
## [23] "PT6"
## [24] "MO_G41"
## [25] "MO_GB1"
## [26] "A1"
## [27] "Liu2012b_Phylo_5"
## [28] "Yamato09_C1"
## [29] "acnaGlo7"
## [30] "MO_G7"
## [31] "Alguacil10_Glo6"
## [32] "Alguacil12b_PARA1"
## [33] "Glo58"
## [34] "lamellosum"
## [35] "VeGlo18"
## [36] "Liu2012a_Ar_1"
## [37] "Div"
## [38] "fennica"
## [39] "Voyriella_parviflora_symbiont_type_1"
## [40] "Para1_OTU2"
## [41] "Whitfield_type_3"
## [42] "Aca"
## [43] "MO_G20"
## [44] "pyriformis"
## [45] "Glo7"
## [46] "MO_GC1"
## [47] "Winther07_B"
## [48] "mosseae"
## [49] "caledonium"
## [50] "brasilianum"
## [51] "Alguacil09b_Glo_G8"
## [52] "Schechter08_Arch1"
## [53] "Wirsel_OTU21"
## [54] "Glo71"
## [55] "spurca"
## [56] "Glo59"
## [57] "Wirsel_OTU16"
## [58] "Glo32"
## [59] "MO_A8"
## [60] "MO_G27"
## [61] "MO_G5"colnames(ssu)## [1] "ssuotu" "MM933M" "MM936M" "MM941M" "MM943M"
## [6] "MM9362CU" "MMS700" "MMS701" "MMS702" "MMS703"
## [11] "MMS704" "MMS705" "MMS706" "MMS707" "MMS708"
## [16] "MMS709" "MMS710" "MMS711" "MMS712" "MMS713"
## [21] "MMS714" "MMS715" "MMS716" "MMS717" "MMS890"
## [26] "MMS891" "MMS900" "MMS901" "MMS902" "MMS903"
## [31] "MMS904" "MMS905" "MMS906" "MMS907" "MMS908"
## [36] "MMS909" "MMS910" "MMS911" "MMS912" "MMS913"
## [41] "MMS914" "MMS915" "MMS916" "MMS917" "MMS918"
## [46] "MMS919" "MMS920" "MMS921" "MMS922" "MMS923"
## [51] "MMS924" "MMS925" "MMS928" "MMS929" "MMS930"
## [56] "MMS931" "MMS932" "MMS938" "MMS939" "MMS942"
## [61] "MMS943" "MMS944" "ssutaxonomy" "ssukingdom" "ssuphylum"
## [66] "ssuclass" "ssuorder" "ssufamily" "ssugenus" "ssuspecies"ssu.save <- ssu
# View(ssu)
# remove Geosiphonaceae
unique(ssu$ssufamily)## [1] "Glomeraceae" NA "Acaulosporaceae"
## [4] "Paraglomeraceae" "Archaeosporaceae" "Ambisporaceae"
## [7] "Claroideoglomeraceae" "Diversisporaceae" "Geosiphonaceae"drop.family <- c("Geosiphonaceae")
ssu <- ssu[-which(ssu$ssufamily %in% drop.family), ]
# remove no blast hits
unique(ssu$ssutaxonomy)## [1] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27"
## [2] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B"
## [3] "No blast hit"
## [4] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp"
## [5] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__Acau16"
## [6] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G47"
## [7] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__laccatum"
## [8] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__acnaGlo2"
## [9] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__MO_Ar1"
## [10] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__sp"
## [11] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__leptoticha"
## [12] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_7"
## [13] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Torrecillas12b_Glo_G5"
## [14] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Sanchez_Castro12b_GLO12"
## [15] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_D"
## [16] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G18"
## [17] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__intraradices"
## [18] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12a_Para_1"
## [19] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Douhan9"
## [20] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Alguacil12b_GLO_G3"
## [21] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo39"
## [22] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo1"
## [23] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Ligrone07_sp"
## [24] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Kuklospora; s__PT6"
## [25] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__sp"
## [26] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G41"
## [27] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__MO_GB1"
## [28] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__A1"
## [29] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Liu2012b_Phylo_5"
## [30] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Yamato09_C1"
## [31] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__acnaGlo7"
## [32] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G7"
## [33] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil10_Glo6"
## [34] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Alguacil12b_PARA1"
## [35] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo58"
## [36] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__lamellosum"
## [37] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__VeGlo18"
## [38] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__Liu2012a_Ar_1"
## [39] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__Div"
## [40] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Ambisporaceae; g__Ambispora; s__fennica"
## [41] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Voyriella_parviflora_symbiont_type_1"
## [42] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__sp"
## [43] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__Para1_OTU2"
## [44] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Whitfield_type_3"
## [45] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Aca"
## [46] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G20"
## [47] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__sp"
## [48] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo7"
## [49] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__MO_GC1"
## [50] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Winther07_B"
## [51] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__mosseae"
## [52] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__caledonium"
## [53] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Paraglomerales; f__Paraglomeraceae; g__Paraglomus; s__brasilianum"
## [54] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Alguacil09b_Glo_G8"
## [55] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Schechter08_Arch1"
## [56] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Archaeosporales; f__Archaeosporaceae; g__Archaeospora; s__Wirsel_OTU21"
## [57] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo71"
## [58] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Diversisporaceae; g__Diversispora; s__spurca"
## [59] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Claroideoglomeraceae; g__Claroideoglomus; s__Glo59"
## [60] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Wirsel_OTU16"
## [61] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__Glo32"
## [62] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Diversisporales; f__Acaulosporaceae; g__Acaulospora; s__MO_A8"
## [63] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G27"
## [64] "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MO_G5"drop <- c("No blast hit")
ssu <- ssu[-which(ssu$ssutaxonomy %in% drop), ]
## alternatively, can use filtering to drop these and avoid rewriting things
ssu <- ssu.save %>% filter(ssufamily != "Geosiphonaceae" & ssutaxonomy != "No blast hit")
ssu <- ssu[complete.cases(ssu), ]
unique(ssu$ssuspecies)## [1] "NES27"
## [2] "MH3_B"
## [3] "sp"
## [4] "Acau16"
## [5] "MO_G47"
## [6] "laccatum"
## [7] "acnaGlo2"
## [8] "MO_Ar1"
## [9] "leptoticha"
## [10] "Whitfield_type_7"
## [11] "Torrecillas12b_Glo_G5"
## [12] "Sanchez_Castro12b_GLO12"
## [13] "Winther07_D"
## [14] "MO_G18"
## [15] "intraradices"
## [16] "Alguacil12a_Para_1"
## [17] "Douhan9"
## [18] "Alguacil12b_GLO_G3"
## [19] "Glo39"
## [20] "Alguacil10_Glo1"
## [21] "Ligrone07_sp"
## [22] "PT6"
## [23] "MO_G41"
## [24] "MO_GB1"
## [25] "A1"
## [26] "Liu2012b_Phylo_5"
## [27] "Yamato09_C1"
## [28] "acnaGlo7"
## [29] "MO_G7"
## [30] "Alguacil10_Glo6"
## [31] "Alguacil12b_PARA1"
## [32] "Glo58"
## [33] "lamellosum"
## [34] "VeGlo18"
## [35] "Liu2012a_Ar_1"
## [36] "Div"
## [37] "fennica"
## [38] "Voyriella_parviflora_symbiont_type_1"
## [39] "Para1_OTU2"
## [40] "Whitfield_type_3"
## [41] "Aca"
## [42] "MO_G20"
## [43] "Glo7"
## [44] "MO_GC1"
## [45] "Winther07_B"
## [46] "mosseae"
## [47] "caledonium"
## [48] "brasilianum"
## [49] "Alguacil09b_Glo_G8"
## [50] "Schechter08_Arch1"
## [51] "Wirsel_OTU21"
## [52] "Glo71"
## [53] "spurca"
## [54] "Glo59"
## [55] "Wirsel_OTU16"
## [56] "Glo32"
## [57] "MO_A8"
## [58] "MO_G27"
## [59] "MO_G5"Add in functional groups
ssul <- melt(ssu, variable.name = "ID", value.name = "ssureads")
str(ssul)## 'data.frame': 17568 obs. of 11 variables:
## $ ssuotu : chr "denovo538" "denovo1813859" "denovo34" "denovo11952" ...
## $ ssutaxonomy: chr "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp" ...
## $ ssukingdom : chr "Fungi" "Fungi" "Fungi" "Fungi" ...
## $ ssuphylum : chr "Glomeromycota" "Glomeromycota" "Glomeromycota" "Glomeromycota" ...
## $ ssuclass : chr "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" ...
## $ ssuorder : chr "Glomerales" "Glomerales" "Glomerales" "Glomerales" ...
## $ ssufamily : chr "Glomeraceae" "Glomeraceae" "Glomeraceae" "Glomeraceae" ...
## $ ssugenus : chr "Glomus" "Glomus" "Glomus" "Glomus" ...
## $ ssuspecies : chr "NES27" "MH3_B" "NES27" "sp" ...
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ssureads : num 5.29 4.33 0 0 3.78 ...OTU richnesses and read abundance for each taxonomic group
Below is in format for graphing
ssugenusREADRICH <- data.frame(ssul %>% group_by(ID, ssugenus) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Read_Abundance_Sample = sum(ssureads)))Make SSU data frame with genus level OTU richness and Read abundance for each genus
ssugenera <- data.frame(ssul %>% group_by(ID, ssugenus) %>% summarise(Genus_OTU_Richness = length(unique(ssuotu[ssureads >
0])), Genus_Read_Abundance = sum(ssureads)))Add metadata into ssu genus
colnames(metassu) # need to edit for appropriate headers## [1] "#SampleID" "BarcodeSequence" "LinkerPrimerSequence"
## [4] "Location" "Project" "Year"
## [7] "Site" "Treat" "Rep"
## [10] "HL" "TSDel" "Description"
## [13] "SampleID2"names(metassu)[names(metassu) == "SampleID2"] <- "ID"
ssugenera <- ssugenera %>% left_join(metassu, by = "ID")
str(ssugenera)## 'data.frame': 488 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssugenus : chr "Acaulospora" "Ambispora" "Archaeospora" "Claroideoglomus" ...
## $ Genus_OTU_Richness : int 2 5 17 15 2 65 2 22 1 5 ...
## $ Genus_Read_Abundance: num 6.66 23.94 157.59 99.62 23.2 ...
## $ #SampleID : int 5 5 5 5 5 5 5 5 3 3 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
## $ LinkerPrimerSequence: logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "E1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 5 3 3 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...Adds meta data to data for box plot
ssugenusREADRICH <- ssugenusREADRICH %>% left_join(metassu, by = "ID")
str(ssugenusREADRICH)## 'data.frame': 488 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssugenus : chr "Acaulospora" "Ambispora" "Archaeospora" "Claroideoglomus" ...
## $ OTU_Richness_Sample : int 2 5 17 15 2 65 2 22 1 5 ...
## $ Read_Abundance_Sample: num 6.66 23.94 157.59 99.62 23.2 ...
## $ #SampleID : int 5 5 5 5 5 5 5 5 3 3 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "E1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 5 3 3 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...filename <- paste0(date, "_ssu_genera_read_rich.csv")
write.csv(ssugenusREADRICH, file = filename, row.names = FALSE)move ID to single line and genus to the length
### For ssufamily, moves ID to single line and genus to the length
LssugenusOTURICH <- dcast(ssugenera, ID ~ ssugenus, value.var = "Genus_OTU_Richness",
fun.aggregate = sum)
LssugenusOTUREAD <- dcast(ssugenera, ID ~ ssugenus, value.var = "Genus_Read_Abundance",
fun.aggregate = sum)
head(LssugenusOTUREAD)## ID Acaulospora Ambispora Archaeospora Claroideoglomus Diversispora
## 1 MM933M 6.6618 23.9417 157.5890 99.6242 23.1980
## 2 MM9362CU 12.2907 33.2473 63.2683 118.7108 18.3924
## 3 MM936M 3.3956 25.9602 67.7783 91.5404 12.2380
## 4 MM941M 21.9872 14.6280 47.8637 93.6373 13.9450
## 5 MM943M 14.2272 23.4654 92.1497 74.7298 22.9824
## 6 MMS700 11.9400 20.1172 90.3423 106.1814 18.6176
## Glomus Kuklospora Paraglomus
## 1 387.3753 5.7632 129.6337
## 2 464.2526 3.1134 172.7572
## 3 406.1616 6.7912 81.1362
## 4 503.8588 11.0910 81.3654
## 5 317.0982 30.6556 103.0184
## 6 429.0388 10.8650 118.3515Genera reads and richness together ###WORKED ON THIS 20180826
REPLACED Functional.Group with ssugenus
Work on changing to species in the next iteration
ssugenusREADRICHlong <- data.frame(ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Acaulospora_Richness = length(unique(ssuotu[ssureads > 0 & ssugenus ==
"Acaulospora"])), Ambispora_Richness = length(unique(ssuotu[ssureads > 0 &
ssugenus == "Ambispora"])), Archaeospora_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Archaeospora"])), Claroideoglomus_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Claroideoglomus"])), Diversispora_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Diversispora"])), Glomus_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Glomus"])), Kuklospora_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Kuklospora"])), Paraglomus_Richness = length(unique(ssuotu[ssureads >
0 & ssugenus == "Paraglomus"])), Acaulospora_Read_Abundance = sum(ssureads[ssugenus ==
"Acaulospora"]), Ambispora_Read_Abundance = sum(ssureads[ssugenus == "Ambispora"]),
Archaeospora_Read_Abundance = sum(ssureads[ssugenus == "Archaeospora"]),
Claroideoglomus_Read_Abundance = sum(ssureads[ssugenus == "Claroideoglomus"]),
Diversispora_Read_Abundance = sum(ssureads[ssugenus == "Diversispora"]),
Glomus_Read_Abundance = sum(ssureads[ssugenus == "Glomus"]), Kuklospora_Read_Abundance = sum(ssureads[ssugenus ==
"Kuklospora"]), Paraglomus_Read_Abundance = sum(ssureads[ssugenus ==
"Paraglomus"]), Read_Abundance_Sample = sum(ssureads)))
## View(ssugenusREADRICHlong)
ssugenusREADRICHlong <- merge(ssugenusREADRICHlong, metassu, by = "ID")
filename <- paste0(date, "_ssu_genus_read_rich.csv")
write.csv(ssugenusREADRICHlong, file = filename, row.names = FALSE)
str(ssugenusREADRICHlong)## 'data.frame': 61 obs. of 31 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ OTU_Richness_Sample : int 130 167 105 121 99 151 158 161 138 148 ...
## $ Acaulospora_Richness : int 2 3 1 4 3 2 4 2 2 3 ...
## $ Ambispora_Richness : int 5 7 5 3 3 5 5 7 8 6 ...
## $ Archaeospora_Richness : int 17 16 12 10 11 15 13 13 8 13 ...
## $ Claroideoglomus_Richness : int 15 20 13 14 11 19 18 21 20 19 ...
## $ Diversispora_Richness : int 2 3 1 1 3 3 2 4 2 1 ...
## $ Glomus_Richness : int 65 85 59 74 51 83 85 80 79 77 ...
## $ Kuklospora_Richness : int 2 2 2 2 3 1 0 2 1 1 ...
## $ Paraglomus_Richness : int 22 31 12 13 14 23 31 32 18 28 ...
## $ Acaulospora_Read_Abundance : num 6.66 12.29 3.4 21.99 14.23 ...
## $ Ambispora_Read_Abundance : num 23.9 33.2 26 14.6 23.5 ...
## $ Archaeospora_Read_Abundance : num 157.6 63.3 67.8 47.9 92.1 ...
## $ Claroideoglomus_Read_Abundance: num 99.6 118.7 91.5 93.6 74.7 ...
## $ Diversispora_Read_Abundance : num 23.2 18.4 12.2 13.9 23 ...
## $ Glomus_Read_Abundance : num 387 464 406 504 317 ...
## $ Kuklospora_Read_Abundance : num 5.76 3.11 6.79 11.09 30.66 ...
## $ Paraglomus_Read_Abundance : num 129.6 172.8 81.1 81.4 103 ...
## $ Read_Abundance_Sample : num 834 886 695 788 678 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ HL : num 0.264 0.739 0.264 0.264 0.739 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...## View(ssugenusREADRICHlong)Non Unifrac distances, try with bray
WORK ON THIS
Permanova
Permanova tests for differences in composition among groups
Reminder - permanova is always based on pairwise distances/dissimilarities.
Can either nest Treat in Year with Year/Treat Or, can restrict permutations within year, using strata=‘Year’ Or, can remove Year, and focus on Treat and other grouping varuables, like Site
Also, dist.j<-vegdist(ssugenusREADRICHlong[,-1], method=‘jaccard’, na.rm=TRUE)
USE BRAY #dist.bc<-vegdist(ssugenusREADRICHlong[,-1], method=‘bray’, na.rm=TRUE)
Part 3 - multivariate community analyses and data visualization
library(dplyr)
library(reshape2)
library(tidyr)
library(vegan)
library(ggplot2) # ggplot resource -http://rpubs.com/collnell/ggplot2
library(tidyverse) # useful packages in 1 - dplyr, ggplot2, tidyr +community matrix at order level with mapping data
otu_map<-read.csv(‘/Users/maltz/Desktop/RdataBreathe/otus_by_class.csv’)
head(ssugenusREADRICHlong)## ID OTU_Richness_Sample Acaulospora_Richness Ambispora_Richness
## 1 MM933M 130 2 5
## 2 MM9362CU 167 3 7
## 3 MM936M 105 1 5
## 4 MM941M 121 4 3
## 5 MM943M 99 3 3
## 6 MMS700 151 2 5
## Archaeospora_Richness Claroideoglomus_Richness Diversispora_Richness
## 1 17 15 2
## 2 16 20 3
## 3 12 13 1
## 4 10 14 1
## 5 11 11 3
## 6 15 19 3
## Glomus_Richness Kuklospora_Richness Paraglomus_Richness
## 1 65 2 22
## 2 85 2 31
## 3 59 2 12
## 4 74 2 13
## 5 51 3 14
## 6 83 1 23
## Acaulospora_Read_Abundance Ambispora_Read_Abundance
## 1 6.6618 23.9417
## 2 12.2907 33.2473
## 3 3.3956 25.9602
## 4 21.9872 14.6280
## 5 14.2272 23.4654
## 6 11.9400 20.1172
## Archaeospora_Read_Abundance Claroideoglomus_Read_Abundance
## 1 157.5890 99.6242
## 2 63.2683 118.7108
## 3 67.7783 91.5404
## 4 47.8637 93.6373
## 5 92.1497 74.7298
## 6 90.3423 106.1814
## Diversispora_Read_Abundance Glomus_Read_Abundance
## 1 23.1980 387.3753
## 2 18.3924 464.2526
## 3 12.2380 406.1616
## 4 13.9450 503.8588
## 5 22.9824 317.0982
## 6 18.6176 429.0388
## Kuklospora_Read_Abundance Paraglomus_Read_Abundance
## 1 5.7632 129.6337
## 2 3.1134 172.7572
## 3 6.7912 81.1362
## 4 11.0910 81.3654
## 5 30.6556 103.0184
## 6 10.8650 118.3515
## Read_Abundance_Sample #SampleID BarcodeSequence LinkerPrimerSequence
## 1 833.7869 5 ATCCCGTATCGATTGG NA
## 2 886.0327 4 GCAACCTTTCGATTGG NA
## 3 695.0015 3 GCAACCTTAGAGTGTG NA
## 4 788.3764 1 GCAACCTTTGGGTGAT NA
## 5 678.3267 2 GCAACCTTTACTGTGC NA
## 6 805.4538 47 GAAGATCCCTCTCAAG NA
## Location Project Year Site Treat Rep HL TSDel Description
## 1 E1 MM-Salvage 2015 WL O 5 0.2644545 <NA> RecipientPre
## 2 D1 MM-Salvage 2015 WL O 4 0.7392792 <NA> RecipientPre
## 3 C1 MM-Salvage 2015 WL O 3 0.2644545 <NA> RecipientPre
## 4 A1 MM-Salvage 2015 WL O 1 0.2644545 <NA> RecipientPre
## 5 B1 MM-Salvage 2015 WL O 2 0.7392792 <NA> RecipientPre
## 6 G6 MM-Salvage 2017 WL D 1 0.4678646 B1A2C RecipientPostcolnames(ssugenusREADRICHlong)## [1] "ID" "OTU_Richness_Sample"
## [3] "Acaulospora_Richness" "Ambispora_Richness"
## [5] "Archaeospora_Richness" "Claroideoglomus_Richness"
## [7] "Diversispora_Richness" "Glomus_Richness"
## [9] "Kuklospora_Richness" "Paraglomus_Richness"
## [11] "Acaulospora_Read_Abundance" "Ambispora_Read_Abundance"
## [13] "Archaeospora_Read_Abundance" "Claroideoglomus_Read_Abundance"
## [15] "Diversispora_Read_Abundance" "Glomus_Read_Abundance"
## [17] "Kuklospora_Read_Abundance" "Paraglomus_Read_Abundance"
## [19] "Read_Abundance_Sample" "#SampleID"
## [21] "BarcodeSequence" "LinkerPrimerSequence"
## [23] "Location" "Project"
## [25] "Year" "Site"
## [27] "Treat" "Rep"
## [29] "HL" "TSDel"
## [31] "Description"Create the comm.grps from ssugenusREADRICHlong
To make the comm.grps I want the first column (col 1), then I don’t want (want to de-select) cols 2 through 22, and then I want the rest of the data = cols 23-30
This didn’t work:
comm.grps<-ssugenusREADRICHlong%>%dplyr::select(‘#SampleID’:Description) #mapping data #heads[1,-c(12:17)] ##This will not work because it uses numeric syntax and the text, only use one or the other (only select; or only numeric -c) comm.grps<-ssugenusREADRICHlong%>%dplyr::select[‘ID’:Description,-c(2:22)] #mapping data
#So I used the longhand way of typing out all the column names that I wanted to keep
`?`(dplyr::select)
comm.grps <- ssugenusREADRICHlong %>% dplyr::select("ID", "Location", "Year",
"Site", "Treat", "Rep", "TSDel", "Description") #mapping data
colnames(comm.grps)## [1] "ID" "Location" "Year" "Site" "Treat"
## [6] "Rep" "TSDel" "Description"WORK ON THIS!!!
Make the comm.mat ##Make the changes for this for species for the next iteration
comm.mat <- ssugenusREADRICHlong %>% dplyr::select(Acaulospora_Richness:Paraglomus_Richness) # community matrix - all but mapping datacomparing ecological communities
diversity vs composition
abundance and richness are univariate response variables used to quantify communities
in multivariate analyses we have these variables for multiple entities
similarly, multivariate analyses have counterparts in univariate stats - t-test, ANOVA, mutliple regression
univariate analyses of diversity
list.files() #shows what is in folder## [1] "2018Sep19_ssu_funcguild_read_rich.csv"
## [2] "2018Sep19_ssu_genera_read_rich.csv"
## [3] "2018Sep19_ssu_genus_read_rich.csv"
## [4] "2018Sep19_ssu_species_read_rich.csv"
## [5] "2018Sep24_ssu_funcguild_read_rich.csv"
## [6] "2018Sep24_ssu_genera_read_rich.csv"
## [7] "2018Sep24_ssu_species_read_rich.csv"
## [8] "2018Sep27_ssu_genera_read_rich.csv"
## [9] "2018Sep27_ssu_genus_read_rich.csv"
## [10] "2018Sep27_ssu_species_read_rich.csv"
## [11] "AMF.Rmd"
## [12] "AMF2.Rmd"
## [13] "AMF4.Rmd"
## [14] "AMF5.1.Rmd"
## [15] "AMFModeling.Rmd"
## [16] "AMFModelingRevise180913.Rmd"
## [17] "AMFModelingSpecies_1.Rmd"
## [18] "AMFModelingTaxa_3.Rmd"
## [19] "AMFModelingTaxa1 2.Rmd"
## [20] "AMFModelingTaxa1.Rmd"
## [21] "AMFModelingTaxa2.Rmd"
## [22] "AMFSubsetRevise180913.Rmd"
## [23] "bray_curtis_nmds_converted.txt"
## [24] "CN.csv"
## [25] "CSS_table_sorted.txt"
## [26] "HL.xlsx"
## [27] "HyphalSlidesDataInput_LD.xlsx"
## [28] "IdentifyingIssues"
## [29] "map_SSU_Salvage_qiimeformat.txt"
## [30] "mapSal.txt"
## [31] "mapSal18924.txt"
## [32] "mapSal18927.txt"
## [33] "rsconnect"
## [34] "SalvageHyphalExtractionSlidesMM.calc_in process.xlsx"
## [35] "SalvageProjectDescription"
## [36] "SSU_species.rtf"
## [37] "Test_cache"
## [38] "Test_files"
## [39] "Test.html"
## [40] "Test.Rmd"
## [41] "Test27_cache"
## [42] "Test27.Rmd"
## [43] "TestSpecies_cache"
## [44] "TestSpecies_files"
## [45] "TestSpecies.html"
## [46] "TestSpecies.Rmd"
## [47] "TestSpeciesCorr2.Rmd"
## [48] "TestSpeciesCorrected_cache"
## [49] "TestSpeciesCorrected_files"
## [50] "TestSpeciesCorrected.html"
## [51] "TestSpeciesCorrected.Rmd"
## [52] "TypesOfAMF"# View(comm.mat)
str(comm.grps)## 'data.frame': 61 obs. of 8 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description: chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...head(comm.grps)## ID Location Year Site Treat Rep TSDel Description
## 1 MM933M E1 2015 WL O 5 <NA> RecipientPre
## 2 MM9362CU D1 2015 WL O 4 <NA> RecipientPre
## 3 MM936M C1 2015 WL O 3 <NA> RecipientPre
## 4 MM941M A1 2015 WL O 1 <NA> RecipientPre
## 5 MM943M B1 2015 WL O 2 <NA> RecipientPre
## 6 MMS700 G6 2017 WL D 1 B1A2C RecipientPoststr(comm.mat)## 'data.frame': 61 obs. of 8 variables:
## $ Acaulospora_Richness : int 2 3 1 4 3 2 4 2 2 3 ...
## $ Ambispora_Richness : int 5 7 5 3 3 5 5 7 8 6 ...
## $ Archaeospora_Richness : int 17 16 12 10 11 15 13 13 8 13 ...
## $ Claroideoglomus_Richness: int 15 20 13 14 11 19 18 21 20 19 ...
## $ Diversispora_Richness : int 2 3 1 1 3 3 2 4 2 1 ...
## $ Glomus_Richness : int 65 85 59 74 51 83 85 80 79 77 ...
## $ Kuklospora_Richness : int 2 2 2 2 3 1 0 2 1 1 ...
## $ Paraglomus_Richness : int 22 31 12 13 14 23 31 32 18 28 ...head(comm.mat)## Acaulospora_Richness Ambispora_Richness Archaeospora_Richness
## 1 2 5 17
## 2 3 7 16
## 3 1 5 12
## 4 4 3 10
## 5 3 3 11
## 6 2 5 15
## Claroideoglomus_Richness Diversispora_Richness Glomus_Richness
## 1 15 2 65
## 2 20 3 85
## 3 13 1 59
## 4 14 1 74
## 5 11 3 51
## 6 19 3 83
## Kuklospora_Richness Paraglomus_Richness
## 1 2 22
## 2 2 31
## 3 2 12
## 4 2 13
## 5 3 14
## 6 1 23does diversity vary across groups?
compute diversity indices
indices <- comm.grps
indices\(richness <- rowSums(comm.mat >0) #indices\)shannon <- diversity(comm.mat, index=‘shannon’)
indices$rarified <- c(rarefy(comm.mat, sample=50)) # rarefied diversity for a given sample size
indices <- comm.grps
indices$richness <- rowSums(comm.mat > 0)
indices$shannon <- diversity(comm.mat, index = "shannon")
View(indices)
# indices$rarified <- c(rarefy(comm.mat, sample=18103)) # rarefied diversity
# for a given sample sizeWORK ON THIS!!!
Depending on what I put first in the model (Treat or Description) The terms are either significant (p-value) or not-significant.
Subset the data
prevpost.mat = comm.mat[comm.grps$Site == "WL", ]
prevpost.grps = comm.grps[comm.grps$Site == "WL", ]
distWL.bc <- vegdist(prevpost.mat, method = "bray", na.rm = TRUE)
WLTreat.div <- adonis2(distWL.bc ~ Description, strata = Treat, data = prevpost.grps,
permutations = 999, method = "bray")
WLTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distWL.bc ~ Description, data = prevpost.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.073503 0.35663 9.9779 0.003 **
## Residual 18 0.132599 0.64337
## Total 19 0.206101 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dvpostWL.mat = comm.mat[(comm.grps$Site == "WL" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostWL.grps = comm.grps[(comm.grps$Site == "WL" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDWL.bc <- vegdist(dvpostWL.mat, method = "bray", na.rm = TRUE)
DWLTreat.div <- adonis2(distDWL.bc ~ Description, strata = Treat, data = dvpostWL.grps,
permutations = 999, method = "bray")
DWLTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDWL.bc ~ Description, data = dvpostWL.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.015150 0.19184 4.0355 0.013 *
## Residual 17 0.063823 0.80816
## Total 18 0.078974 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1prevpostPS.mat = comm.mat[comm.grps$Site == "PS", ]
prevpostPS.grps = comm.grps[comm.grps$Site == "PS", ]
distPS.bc <- vegdist(prevpostPS.mat, method = "bray", na.rm = TRUE)
PSTreat.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = prevpostPS.grps,
permutations = 999, method = "bray")
PSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distPS.bc ~ Description, data = prevpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.012597 0.17068 3.4988 0.031 *
## Residual 17 0.061209 0.82932
## Total 18 0.073806 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(prevpostPS.grps)## [1] 19 8View(prevpostPS.grps)
dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps,
permutations = 999, method = "bray")
DPSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.008275 0.1119 2.268 0.079 .
## Residual 18 0.065677 0.8881
## Total 19 0.073953 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(dvpostPS.grps)## [1] 20 8postTreatPS.mat = comm.mat[comm.grps$Site == "PS", ]
postTreatPS.grps = comm.grps[comm.grps$Site == "PS", ]
distTPS.bc <- vegdist(postTreatPS.mat, method = "bray", na.rm = TRUE)
View(postTreatPS.grps)
contr.mineCD <- function(...) cbind(c(-1, 1, 0, 0, 0, 0))
PSTreatCD <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineCD")
PSTreatCD## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCD")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.034883 0.47264 2.3302 0.022 *
## Residual 13 0.038923 0.52736
## Total 18 0.073806 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1contr.mineCT <- function(...) cbind(c(-1, 0, 1, 0, 0, 0))
PSTreatCT <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineCT")
PSTreatCT## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineCT")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.034883 0.47264 2.3302 0.018 *
## Residual 13 0.038923 0.52736
## Total 18 0.073806 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1contr.mineOC <- function(...) cbind(c(-1, 0, 0, 0, 0, 1))
PSTreatOC <- adonis2(distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999,
method = "bray", contr.unordered = "contr.mineOC")
PSTreatOC## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distTPS.bc ~ Treat, data = postTreatPS.grps, permutations = 999, method = "bray", contr.unordered = "contr.mineOC")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.034883 0.47264 2.3302 0.028 *
## Residual 13 0.038923 0.52736
## Total 18 0.073806 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1postTreatPS.div <- adonis2(distPS.bc ~ Description, strata = Treat, data = postTreatPS.grps,
permutations = 999, method = "bray")
postTreatPS.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distPS.bc ~ Description, data = postTreatPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.012597 0.17068 3.4988 0.025 *
## Residual 17 0.061209 0.82932
## Total 18 0.073806 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(postTreatPS.grps)## [1] 19 8View(postTreatPS.grps)
dvpostPS.mat = comm.mat[(comm.grps$Site == "PS" & comm.grps$Description == "RecipientPost") |
comm.grps$Site == "DS", ]
dvpostPS.grps = comm.grps[(comm.grps$Site == "PS" & comm.grps$Description ==
"RecipientPost") | comm.grps$Site == "DS", ]
distDPS.bc <- vegdist(dvpostPS.mat, method = "bray", na.rm = TRUE)
DPSTreat.div <- adonis2(distDPS.bc ~ Description, strata = Treat, data = dvpostPS.grps,
permutations = 999, method = "bray")
DPSTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = distDPS.bc ~ Description, data = dvpostPS.grps, permutations = 999, method = "bray", strata = Treat)
## Df SumOfSqs R2 F Pr(>F)
## Description 1 0.008275 0.1119 2.268 0.09 .
## Residual 18 0.065677 0.8881
## Total 19 0.073953 1.0000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1dim(dvpostPS.grps)## [1] 20 8set.seed(304)
`?`(vegdist)
## are the results the same with other (non evolutionary) dissimiarlity
## indices? dist.j<-vegdist(ssugenusREADRICHlong[,-1], method='jaccard')
## dist.bc<-vegdist(ssugenusREADRICHlong[,-1], method='bray', na.rm=TRUE)
## dist.j<-vegdist(comm.mat, method='jaccard')
dist.bc <- vegdist(comm.mat, method = "bray", na.rm = TRUE)
AMTreat.div <- adonis2(dist.bc ~ Treat, data = comm.grps, permutations = 999,
method = "bray")
AMTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.09754 0.24065 3.486 0.001 ***
## Residual 55 0.30779 0.75935
## Total 60 0.40533 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AMTD.div <- adonis2(dist.bc ~ Treat + Description, data = comm.grps, permutations = 999,
method = "bray")
AMTD.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat + Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.09754 0.24065 3.5118 0.003 **
## Description 1 0.00782 0.01928 1.4071 0.227
## Residual 54 0.29997 0.74007
## Total 60 0.40533 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1AMTSD.div <- adonis2(dist.bc ~ Treat + Site * Description, data = comm.grps,
permutations = 999, method = "bray")
AMTSD.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat + Site * Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.09754 0.24065 3.9234 0.001 ***
## Site 3 0.03184 0.07857 2.1348 0.057 .
## Site:Description 2 0.02733 0.06742 2.7480 0.052 .
## Residual 50 0.24861 0.61337
## Total 60 0.40533 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# AMD_ST_J.div<-adonis2(dist.bc~Description+Site*Treat, data=comm.grps,
# permutations = 999, method='bray') AMD_ST_J.div
AMD_ST_BC.div <- adonis2(dist.bc ~ Description + Site + Treat, data = comm.grps,
permutations = 999, method = "bray")
AMD_ST_BC.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Description + Site + Treat, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Description 2 0.07616 0.18790 7.1762 0.002 **
## Site 2 0.02444 0.06030 2.3030 0.065 .
## Treat 4 0.02878 0.07101 1.3559 0.230
## Residual 52 0.27594 0.68079
## Total 60 0.40533 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Depending on what comes first in the model, we see different significance of these factors
AMTreat.div <- adonis2(dist.bc ~ Treat * Description, data = comm.grps, permutations = 999,
method = "bray")
AMTreat.div## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dist.bc ~ Treat * Description, data = comm.grps, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Treat 5 0.09754 0.24065 3.5118 0.002 **
## Description 1 0.00782 0.01928 1.4071 0.231
## Residual 54 0.29997 0.74007
## Total 60 0.40533 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1amfMDS<-metaMDS(comm.mat, distance="bray", k=2, trymax=35, autotransform=TRUE) ##k is the number of dimensions## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1032977
## Run 1 stress 0.1190699
## Run 2 stress 0.1180244
## Run 3 stress 0.129419
## Run 4 stress 0.1026719
## ... New best solution
## ... Procrustes: rmse 0.02425996 max resid 0.1500037
## Run 5 stress 0.1144604
## Run 6 stress 0.1221095
## Run 7 stress 0.1212653
## Run 8 stress 0.1283063
## Run 9 stress 0.1080377
## Run 10 stress 0.1182926
## Run 11 stress 0.1220324
## Run 12 stress 0.1245569
## Run 13 stress 0.1218783
## Run 14 stress 0.1067808
## Run 15 stress 0.1033579
## Run 16 stress 0.1062276
## Run 17 stress 0.1033075
## Run 18 stress 0.1358752
## Run 19 stress 0.1033716
## Run 20 stress 0.1257308
## Run 21 stress 0.1144457
## Run 22 stress 0.1040688
## Run 23 stress 0.1342131
## Run 24 stress 0.1206015
## Run 25 stress 0.1087087
## Run 26 stress 0.1164982
## Run 27 stress 0.4024507
## Run 28 stress 0.1224657
## Run 29 stress 0.1147626
## Run 30 stress 0.1223208
## Run 31 stress 0.1304137
## Run 32 stress 0.1166109
## Run 33 stress 0.1023312
## ... New best solution
## ... Procrustes: rmse 0.01416885 max resid 0.06962781
## Run 34 stress 0.1036938
## Run 35 stress 0.120806
## *** No convergence -- monoMDS stopping criteria:
## 35: stress ratio > sratmaxamfMDS ##metaMDS takes eaither a distance matrix or your community matrix (then requires method for 'distance=')##
## Call:
## metaMDS(comm = comm.mat, distance = "bray", k = 2, trymax = 35, autotransform = TRUE)
##
## global Multidimensional Scaling using monoMDS
##
## Data: wisconsin(sqrt(comm.mat))
## Distance: bray
##
## Dimensions: 2
## Stress: 0.1023312
## Stress type 1, weak ties
## No convergent solutions - best solution after 35 tries
## Scaling: centring, PC rotation, halfchange scaling
## Species: expanded scores based on 'wisconsin(sqrt(comm.mat))'stressplot(amfMDS)
#install.packages('ggplot2') ##plotting package
library(ggplot2)
NMDS1 <- amfMDS$points[,1] ##also found using scores(amfMDS)
NMDS2 <- amfMDS$points[,2]
?cbind
amf.plot<-cbind(comm.grps, NMDS1, NMDS2, comm.mat)
p<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Year))+
geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points
stat_ellipse(type='t',size =1)+ ##draws 95% confidence interval ellipses
theme_minimal()
p
plot<-ggplot(amf.plot, aes(NMDS1, NMDS2, color=Treat))+
stat_ellipse(type='t',size =1)+
theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+
annotate("text", x=min(NMDS1), y=min(NMDS2), label=paste('Stress =',round(amfMDS$stress,3))) #add stress to plot
plot
Subset the community data matrix Number of samples that are left
summary(comm.grps)## ID Location Year Site
## MM933M : 1 Length:61 Length:61 Length:61
## MM936M : 1 Class :character Class :character Class :character
## MM941M : 1 Mode :character Mode :character Mode :character
## MM943M : 1
## MM9362CU: 1
## MMS700 : 1
## (Other) :55
## Treat Rep TSDel Description
## C: 9 1:19 Length:61 Length:61
## D: 9 2:19 Class :character Class :character
## F: 9 3:17 Mode :character Mode :character
## O:17 4: 3
## S: 8 5: 2
## T: 9 6: 1
## class(comm.mat)## [1] "data.frame"# subset the comm.mat and comm.grps x=comm.mat[comm.grps$Year == '2015',]
# y=comm.grps[comm.grps$Year == '2015',]commented out some difficult lines
# what is x?
# SSMDS<-metaMDS(x, distance='bray', k = 2, trymax = 35, autotransform =
# TRUE) ##k is the number of dimensions SSMDS ##metaMDS takes eaither a
# distance matrix or your community matrix (then requires method for
# 'distance=')
# stressplot(SSMDS)
# install.packages('ggplot2') ##plotting package
library(ggplot2)
# NMDS1 <- SSMDS$points[,1] ##also found using scores(amfMDS) NMDS2 <-
# SSMDS$points[,2]
`?`(cbind)The code was stopping at this point, and I created a new chunk below
SSamf.plot<-cbind(y, NMDS1, NMDS2, x)
Sp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=7)+##separates overlapping points # stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() Sp
SSp<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Site))+ geom_point(position=position_jitter(.1), shape=3)+##separates overlapping points stat_ellipse(type=‘t’,size =1)+ ##draws 95% confidence interval ellipses theme_minimal() SSp
plot<-ggplot(SSamf.plot, aes(NMDS1, NMDS2, color=Treat))+ stat_ellipse(type=‘t’,size =1)+ theme_minimal()+geom_text(data=amf.plot,aes(NMDS1, NMDS2, label=Site), position=position_jitter(.35))+ annotate(“text”, x=min(NMDS1), y=min(NMDS2), label=paste(‘Stress =’,round(amfMDS$stress,3))) #add stress to plot plot
ad.bc<-adonis2(dist.bc~Treat+Metadata1+Metadata, data=grps, permutations=1000)
ad.bc
Examining the statistic support for increased rhizophilic richness in recipientpost plots
ssu.fg <- ssu %>% mutate(Functional.Group = case_when(ssufamily %in% c("Glomeraceae",
"Claroideoglomeraceae", "Paraglomeraceae") ~ "Rhizophilic", ssufamily %in%
c("Gigasporaceae", "Diversisporaceae") ~ "Edaphophilic", ssufamily %in%
c("Ambisporaceae", "Archaeosporaceae", "Acaulosporaceae") ~ "Ancestral"))
ssul <- melt(ssu.fg, variable.name = "ID", value.name = "ssureads")
str(ssul)## 'data.frame': 17568 obs. of 12 variables:
## $ ssuotu : chr "denovo538" "denovo1813859" "denovo34" "denovo11952" ...
## $ ssutaxonomy : chr "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__MH3_B" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__NES27" "k__Fungi; p__Glomeromycota; c__Glomeromycetes; o__Glomerales; f__Glomeraceae; g__Glomus; s__sp" ...
## $ ssukingdom : chr "Fungi" "Fungi" "Fungi" "Fungi" ...
## $ ssuphylum : chr "Glomeromycota" "Glomeromycota" "Glomeromycota" "Glomeromycota" ...
## $ ssuclass : chr "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" "Glomeromycetes" ...
## $ ssuorder : chr "Glomerales" "Glomerales" "Glomerales" "Glomerales" ...
## $ ssufamily : chr "Glomeraceae" "Glomeraceae" "Glomeraceae" "Glomeraceae" ...
## $ ssugenus : chr "Glomus" "Glomus" "Glomus" "Glomus" ...
## $ ssuspecies : chr "NES27" "MH3_B" "NES27" "sp" ...
## $ Functional.Group: chr "Rhizophilic" "Rhizophilic" "Rhizophilic" "Rhizophilic" ...
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ ssureads : num 5.29 4.33 0 0 3.78 ...ssuguildREADRICH <- data.frame(ssul %>% group_by(ID, Functional.Group) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Read_Abundance_Sample = sum(ssureads)))
ssufamily <- data.frame(ssul %>% group_by(ID, ssufamily) %>% summarise(Family_OTU_Richness = length(unique(ssuotu[ssureads >
0])), Family_Read_Abundance = sum(ssureads)))
colnames(metassu) # need to edit for appropriate headers## [1] "#SampleID" "BarcodeSequence" "LinkerPrimerSequence"
## [4] "Location" "Project" "Year"
## [7] "Site" "Treat" "Rep"
## [10] "HL" "TSDel" "Description"
## [13] "ID"names(metassu)[names(metassu) == "SampleID2"] <- "ID"
ssufamily <- ssufamily %>% left_join(metassu, by = "ID")
str(ssufamily)## 'data.frame': 427 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM933M" ...
## $ ssufamily : chr "Acaulosporaceae" "Ambisporaceae" "Archaeosporaceae" "Claroideoglomeraceae" ...
## $ Family_OTU_Richness : int 4 5 17 15 2 65 22 3 5 12 ...
## $ Family_Read_Abundance: num 12.4 23.9 157.6 99.6 23.2 ...
## $ #SampleID : int 5 5 5 5 5 5 5 3 3 3 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "E1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 5 5 5 5 3 3 3 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...ssuguildREADRICH <- ssuguildREADRICH %>% left_join(metassu, by = "ID")
str(ssuguildREADRICH)## 'data.frame': 183 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICH, file = filename, row.names = FALSE)
LssufamilyOTURICH <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_OTU_Richness",
fun.aggregate = sum)
LssufamilyOTUREAD <- dcast(ssufamily, ID ~ ssufamily, value.var = "Family_Read_Abundance",
fun.aggregate = sum)
head(LssufamilyOTUREAD)## ID Acaulosporaceae Ambisporaceae Archaeosporaceae
## 1 MM933M 12.4250 23.9417 157.5890
## 2 MM9362CU 15.4041 33.2473 63.2683
## 3 MM936M 10.1868 25.9602 67.7783
## 4 MM941M 33.0782 14.6280 47.8637
## 5 MM943M 44.8828 23.4654 92.1497
## 6 MMS700 22.8050 20.1172 90.3423
## Claroideoglomeraceae Diversisporaceae Glomeraceae Paraglomeraceae
## 1 99.6242 23.1980 387.3753 129.6337
## 2 118.7108 18.3924 464.2526 172.7572
## 3 91.5404 12.2380 406.1616 81.1362
## 4 93.6373 13.9450 503.8588 81.3654
## 5 74.7298 22.9824 317.0982 103.0184
## 6 106.1814 18.6176 429.0388 118.3515ssuguildREADRICHlong <- data.frame(ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group ==
"Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Rhizophilic"])), Ancestral_Read_Abundance = sum(ssureads[Functional.Group ==
"Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Rhizophilic"]), Read_Abundance_Sample = sum(ssureads)))
View(ssuguildREADRICHlong)
ssuguildREADRICHlong <- merge(ssuguildREADRICHlong, metassu, by = "ID")
filename <- paste0(date, "_ssu_funcguild_read_rich.csv")
write.csv(ssuguildREADRICHlong, file = filename, row.names = FALSE)
str(ssuguildREADRICHlong)## 'data.frame': 61 obs. of 21 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ OTU_Richness_Sample : int 130 167 105 121 99 151 158 161 138 148 ...
## $ Ancestral_Richness : int 26 28 20 19 20 23 22 24 19 23 ...
## $ Edaphophilic_Richness : int 2 3 1 1 3 3 2 4 2 1 ...
## $ Rhizophilic_Richness : int 102 136 84 101 76 125 134 133 117 124 ...
## $ Ancestral_Read_Abundance : num 194 111.9 103.9 95.6 160.5 ...
## $ Edaphophilic_Read_Abundance: num 23.2 18.4 12.2 13.9 23 ...
## $ Rhizophilic_Read_Abundance : num 617 756 579 679 495 ...
## $ Read_Abundance_Sample : num 834 886 695 788 678 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ HL : num 0.264 0.739 0.264 0.264 0.739 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...ssuguildREAD <- ssul %>% group_by(ID) %>% summarise(Ancestral_Read_Abundance = sum(ssureads[Functional.Group ==
"Ancestral"]), Edaphophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Edaphophilic"]), Rhizophilic_Read_Abundance = sum(ssureads[Functional.Group ==
"Rhizophilic"]))
ssuguildRICH <- ssul %>% group_by(ID) %>% summarise(OTU_Richness_Sample = length(unique(ssuotu[ssureads >
0])), Ancestral_Richness = length(unique(ssuotu[ssureads > 0 & Functional.Group ==
"Ancestral"])), Edaphophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Edaphophilic"])), Rhizophilic_Richness = length(unique(ssuotu[ssureads >
0 & Functional.Group == "Rhizophilic"])))
head(ssuguildRICH)## # A tibble: 6 x 5
## ID OTU_Richness_Sample Ancestral_Richness Edaphophilic_Richness
## <fct> <int> <int> <int>
## 1 MM933M 130 26 2
## 2 MM936M 105 20 1
## 3 MM941M 121 19 1
## 4 MM943M 99 20 3
## 5 MM9362CU 167 28 3
## 6 MMS700 151 23 3
## # ... with 1 more variable: Rhizophilic_Richness <int>ssuguildRICH <- merge(ssuguildRICH, metassu, by = "ID")
ssuguildREAD <- merge(ssuguildREAD, metassu, by = "ID")
str(ssuguildRICH)## 'data.frame': 61 obs. of 17 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ OTU_Richness_Sample : int 130 167 105 121 99 151 158 161 138 148 ...
## $ Ancestral_Richness : int 26 28 20 19 20 23 22 24 19 23 ...
## $ Edaphophilic_Richness: int 2 3 1 1 3 3 2 4 2 1 ...
## $ Rhizophilic_Richness : int 102 136 84 101 76 125 134 133 117 124 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ HL : num 0.264 0.739 0.264 0.264 0.739 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...str(ssuguildREAD)## 'data.frame': 61 obs. of 16 variables:
## $ ID : Factor w/ 61 levels "MM933M","MM936M",..: 1 5 2 3 4 6 7 8 9 10 ...
## $ Ancestral_Read_Abundance : num 194 111.9 103.9 95.6 160.5 ...
## $ Edaphophilic_Read_Abundance: num 23.2 18.4 12.2 13.9 23 ...
## $ Rhizophilic_Read_Abundance : num 617 756 579 679 495 ...
## $ #SampleID : int 5 4 3 1 2 47 48 49 50 51 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "GCAACCTTTCGATTGG" "GCAACCTTAGAGTGTG" "GCAACCTTTGGGTGAT" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "D1" "C1" "A1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 2 2 2 2 2 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 4 3 1 2 1 2 3 1 2 ...
## $ HL : num 0.264 0.739 0.264 0.264 0.739 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...y <- ssuguildRICH$Rhizophilic_Richness
str(y)## int [1:61] 102 136 84 101 76 125 134 133 117 124 ...D <- ssuguildRICH$Description
S <- ssuguildRICH$Site
length(S)## [1] 61length(D)## [1] 61`?`(lm)
# lm(y~D)
SRhizRich <- lm(y ~ S)
str(SRhizRich)## List of 13
## $ coefficients : Named num [1:4] 125.4 12.6 6.18 3.2
## ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ residuals : Named num [1:61] -26.6 7.4 -44.6 -27.6 -52.6 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 29.7 11.9 -6.4 -44.7 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 4
## $ fitted.values: Named num [1:61] 129 129 129 129 129 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:4] 0 1 1 1
## $ qr :List of 5
## ..$ qr : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## .. ..- attr(*, "assign")= int [1:4] 0 1 1 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ S: chr "contr.treatment"
## ..$ qraux: num [1:4] 1.13 1.07 1.11 1.07
## ..$ pivot: int [1:4] 1 2 3 4
## ..$ tol : num 1e-07
## ..$ rank : int 4
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 57
## $ contrasts :List of 1
## ..$ S: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## $ call : language lm(formula = y ~ S)
## $ terms :Classes 'terms', 'formula' language y ~ S
## .. ..- attr(*, "variables")= language list(y, S)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. ..$ : chr "S"
## .. ..- attr(*, "term.labels")= chr "S"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S
## .. .. ..- attr(*, "variables")= language list(y, S)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. .. ..$ : chr "S"
## .. .. ..- attr(*, "term.labels")= chr "S"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## - attr(*, "class")= chr "lm"DRhizRich <- lm(y ~ D)
str(DRhizRich)## List of 13
## $ coefficients : Named num [1:3] 125.4 11.44 -8.98
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ residuals : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
## $ rank : int 3
## $ fitted.values: Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:3] 0 1 1
## $ qr :List of 5
## ..$ qr : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## .. ..- attr(*, "assign")= int [1:3] 0 1 1
## .. ..- attr(*, "contrasts")=List of 1
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:3] 1.13 1.18 1.12
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-07
## ..$ rank : int 3
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 58
## $ contrasts :List of 1
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ D)
## $ terms :Classes 'terms', 'formula' language y ~ D
## .. ..- attr(*, "variables")= language list(y, D)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. ..$ : chr "D"
## .. ..- attr(*, "term.labels")= chr "D"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ D
## .. .. ..- attr(*, "variables")= language list(y, D)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. .. ..$ : chr "D"
## .. .. ..- attr(*, "term.labels")= chr "D"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## - attr(*, "class")= chr "lm"DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)## List of 13
## $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.1 ...
## ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ residuals : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 29.7 11.9 -6.4 58.6 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 7
## $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## $ qr :List of 5
## ..$ qr : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## .. ..- attr(*, "contrasts")=List of 2
## .. .. ..$ S: chr "contr.treatment"
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:12] 1.13 1.07 1.11 1.07 1.16 ...
## ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
## ..$ tol : num 1e-07
## ..$ rank : int 7
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 54
## $ contrasts :List of 2
## ..$ S: chr "contr.treatment"
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 2
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ S * D)
## $ terms :Classes 'terms', 'formula' language y ~ S * D
## .. ..- attr(*, "variables")= language list(y, S, D)
## .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## $ model :'data.frame': 61 obs. of 3 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S * D
## .. .. ..- attr(*, "variables")= language list(y, S, D)
## .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## - attr(*, "class")= chr "lm"summary(DSRhizRich)##
## Call:
## lm(formula = y ~ S * D)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.533 -7.533 2.000 7.071 35.167
##
## Coefficients: (5 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 6.941 18.066 < 2e-16 ***
## SHH 8.600 12.986 0.662 0.510616
## SPS -3.400 10.412 -0.327 0.745266
## SWL -18.567 9.398 -1.976 0.053328 .
## DRecipientPost 31.095 7.573 4.106 0.000137 ***
## DRecipientPre NA NA NA NA
## SHH:DRecipientPost -26.562 13.924 -1.908 0.061754 .
## SPS:DRecipientPost -18.962 11.560 -1.640 0.106764
## SWL:DRecipientPost NA NA NA NA
## SHH:DRecipientPre NA NA NA NA
## SPS:DRecipientPre NA NA NA NA
## SWL:DRecipientPre NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared: 0.3019, Adjusted R-squared: 0.2243
## F-statistic: 3.892 on 6 and 54 DF, p-value: 0.002662y <- ssuguildRICH$Rhizophilic_Richness
str(y)## int [1:61] 102 136 84 101 76 125 134 133 117 124 ...D <- ssuguildRICH$Description
S <- ssuguildRICH$Site
length(S)## [1] 61length(D)## [1] 61`?`(lm)
glm(y ~ D)##
## Call: glm(formula = y ~ D)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.400 11.441 -8.983
##
## Degrees of Freedom: 60 Total (i.e. Null); 58 Residual
## Null Deviance: 18630
## Residual Deviance: 14470 AIC: 514.7G_SRhizRich <- glm(y ~ S)
str(G_SRhizRich)## List of 30
## $ coefficients : Named num [1:4] 125.4 12.6 6.18 3.2
## ..- attr(*, "names")= chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ residuals : Named num [1:61] -26.6 7.4 -44.6 -27.6 -52.6 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ fitted.values : Named num [1:61] 129 129 129 129 129 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 29.7 11.9 -6.4 -44.7 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ R : num [1:4, 1:4] -7.81 0 0 0 -2.18 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## $ rank : int 4
## $ qr :List of 5
## ..$ qr : num [1:61, 1:4] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:4] "(Intercept)" "SHH" "SPS" "SWL"
## ..$ rank : int 4
## ..$ qraux: num [1:4] 1.13 1.07 1.11 1.07
## ..$ pivot: int [1:4] 1 2 3 4
## ..$ tol : num 1e-11
## ..- attr(*, "class")= chr "qr"
## $ family :List of 11
## ..$ family : chr "gaussian"
## ..$ link : chr "identity"
## ..$ linkfun :function (mu)
## ..$ linkinv :function (eta)
## ..$ variance :function (mu)
## ..$ dev.resids:function (y, mu, wt)
## ..$ aic :function (y, n, mu, wt, dev)
## ..$ mu.eta :function (eta)
## ..$ initialize: expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$lin| __truncated__
## ..$ validmu :function (mu)
## ..$ valideta :function (eta)
## ..- attr(*, "class")= chr "family"
## $ linear.predictors: Named num [1:61] 129 129 129 129 129 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ deviance : num 17571
## $ aic : num 529
## $ null.deviance : num 18634
## $ iter : int 2
## $ weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ prior.weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ df.residual : int 57
## $ df.null : int 60
## $ y : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ converged : logi TRUE
## $ boundary : logi FALSE
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S
## .. .. ..- attr(*, "variables")= language list(y, S)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. .. ..$ : chr "S"
## .. .. ..- attr(*, "term.labels")= chr "S"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ call : language glm(formula = y ~ S)
## $ formula :Class 'formula' language y ~ S
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## $ terms :Classes 'terms', 'formula' language y ~ S
## .. ..- attr(*, "variables")= language list(y, S)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "S"
## .. .. .. ..$ : chr "S"
## .. ..- attr(*, "term.labels")= chr "S"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "S"
## $ data :<environment: R_GlobalEnv>
## $ offset : NULL
## $ control :List of 3
## ..$ epsilon: num 1e-08
## ..$ maxit : num 25
## ..$ trace : logi FALSE
## $ method : chr "glm.fit"
## $ contrasts :List of 1
## ..$ S: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## - attr(*, "class")= chr [1:2] "glm" "lm"G_DRhizRich <- glm(y ~ D)
str(G_DRhizRich)## List of 30
## $ coefficients : Named num [1:3] 125.4 11.44 -8.98
## ..- attr(*, "names")= chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ residuals : Named num [1:61] -14.4 19.6 -32.4 -15.4 -40.4 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ fitted.values : Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 62.3 16.9 -13.5 -38.5 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "DRecipientPost" "DRecipientPre" "" ...
## $ R : num [1:3, 1:3] -7.81 0 0 -5.63 3.5 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## $ rank : int 3
## $ qr :List of 5
## ..$ qr : num [1:61, 1:3] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:3] "(Intercept)" "DRecipientPost" "DRecipientPre"
## ..$ rank : int 3
## ..$ qraux: num [1:3] 1.13 1.18 1.12
## ..$ pivot: int [1:3] 1 2 3
## ..$ tol : num 1e-11
## ..- attr(*, "class")= chr "qr"
## $ family :List of 11
## ..$ family : chr "gaussian"
## ..$ link : chr "identity"
## ..$ linkfun :function (mu)
## ..$ linkinv :function (eta)
## ..$ variance :function (mu)
## ..$ dev.resids:function (y, mu, wt)
## ..$ aic :function (y, n, mu, wt, dev)
## ..$ mu.eta :function (eta)
## ..$ initialize: expression({ n <- rep.int(1, nobs) if (is.null(etastart) && is.null(start) && is.null(mustart) && ((family$lin| __truncated__
## ..$ validmu :function (mu)
## ..$ valideta :function (eta)
## ..- attr(*, "class")= chr "family"
## $ linear.predictors: Named num [1:61] 116 116 116 116 116 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ deviance : num 14472
## $ aic : num 515
## $ null.deviance : num 18634
## $ iter : int 2
## $ weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ prior.weights : Named num [1:61] 1 1 1 1 1 1 1 1 1 1 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ df.residual : int 58
## $ df.null : int 60
## $ y : Named int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ converged : logi TRUE
## $ boundary : logi FALSE
## $ model :'data.frame': 61 obs. of 2 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ D
## .. .. ..- attr(*, "variables")= language list(y, D)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. .. ..$ : chr "D"
## .. .. ..- attr(*, "term.labels")= chr "D"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ call : language glm(formula = y ~ D)
## $ formula :Class 'formula' language y ~ D
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## $ terms :Classes 'terms', 'formula' language y ~ D
## .. ..- attr(*, "variables")= language list(y, D)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "y" "D"
## .. .. .. ..$ : chr "D"
## .. ..- attr(*, "term.labels")= chr "D"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "character"
## .. .. ..- attr(*, "names")= chr [1:2] "y" "D"
## $ data :<environment: R_GlobalEnv>
## $ offset : NULL
## $ control :List of 3
## ..$ epsilon: num 1e-08
## ..$ maxit : num 25
## ..$ trace : logi FALSE
## $ method : chr "glm.fit"
## $ contrasts :List of 1
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 1
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## - attr(*, "class")= chr [1:2] "glm" "lm"summary(G_DRhizRich)##
## Call:
## glm(formula = y ~ D)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -40.841 -8.400 3.159 10.159 26.159
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 7.064 17.751 <2e-16 ***
## DRecipientPost 11.441 7.455 1.535 0.13
## DRecipientPre -8.983 8.408 -1.068 0.29
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 249.5173)
##
## Null deviance: 18634 on 60 degrees of freedom
## Residual deviance: 14472 on 58 degrees of freedom
## AIC: 514.73
##
## Number of Fisher Scoring iterations: 2DSRhizRich <- lm(y ~ S * D)
str(DSRhizRich)## List of 13
## $ coefficients : Named num [1:12] 125.4 8.6 -3.4 -18.6 31.1 ...
## ..- attr(*, "names")= chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ residuals : Named num [1:61] -4.83 29.17 -22.83 -5.83 -30.83 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ effects : Named num [1:61] -1030.1 29.7 11.9 -6.4 58.6 ...
## ..- attr(*, "names")= chr [1:61] "(Intercept)" "SHH" "SPS" "SWL" ...
## $ rank : int 7
## $ fitted.values: Named num [1:61] 107 107 107 107 107 ...
## ..- attr(*, "names")= chr [1:61] "1" "2" "3" "4" ...
## $ assign : int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## $ qr :List of 5
## ..$ qr : num [1:61, 1:12] -7.81 0.128 0.128 0.128 0.128 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:61] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:12] "(Intercept)" "SHH" "SPS" "SWL" ...
## .. ..- attr(*, "assign")= int [1:12] 0 1 1 1 2 2 3 3 3 3 ...
## .. ..- attr(*, "contrasts")=List of 2
## .. .. ..$ S: chr "contr.treatment"
## .. .. ..$ D: chr "contr.treatment"
## ..$ qraux: num [1:12] 1.13 1.07 1.11 1.07 1.16 ...
## ..$ pivot: int [1:12] 1 2 3 4 5 7 8 6 9 10 ...
## ..$ tol : num 1e-07
## ..$ rank : int 7
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 54
## $ contrasts :List of 2
## ..$ S: chr "contr.treatment"
## ..$ D: chr "contr.treatment"
## $ xlevels :List of 2
## ..$ S: chr [1:4] "DS" "HH" "PS" "WL"
## ..$ D: chr [1:3] "Donor" "RecipientPost" "RecipientPre"
## $ call : language lm(formula = y ~ S * D)
## $ terms :Classes 'terms', 'formula' language y ~ S * D
## .. ..- attr(*, "variables")= language list(y, S, D)
## .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(y, S, D)
## .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## $ model :'data.frame': 61 obs. of 3 variables:
## ..$ y: int [1:61] 102 136 84 101 76 125 134 133 117 124 ...
## ..$ S: chr [1:61] "WL" "WL" "WL" "WL" ...
## ..$ D: chr [1:61] "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language y ~ S * D
## .. .. ..- attr(*, "variables")= language list(y, S, D)
## .. .. ..- attr(*, "factors")= int [1:3, 1:3] 0 1 0 0 0 1 0 1 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:3] "y" "S" "D"
## .. .. .. .. ..$ : chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "term.labels")= chr [1:3] "S" "D" "S:D"
## .. .. ..- attr(*, "order")= int [1:3] 1 1 2
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(y, S, D)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:3] "numeric" "character" "character"
## .. .. .. ..- attr(*, "names")= chr [1:3] "y" "S" "D"
## - attr(*, "class")= chr "lm"summary(DSRhizRich)##
## Call:
## lm(formula = y ~ S * D)
##
## Residuals:
## Min 1Q Median 3Q Max
## -42.533 -7.533 2.000 7.071 35.167
##
## Coefficients: (5 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 125.400 6.941 18.066 < 2e-16 ***
## SHH 8.600 12.986 0.662 0.510616
## SPS -3.400 10.412 -0.327 0.745266
## SWL -18.567 9.398 -1.976 0.053328 .
## DRecipientPost 31.095 7.573 4.106 0.000137 ***
## DRecipientPre NA NA NA NA
## SHH:DRecipientPost -26.562 13.924 -1.908 0.061754 .
## SPS:DRecipientPost -18.962 11.560 -1.640 0.106764
## SWL:DRecipientPost NA NA NA NA
## SHH:DRecipientPre NA NA NA NA
## SPS:DRecipientPre NA NA NA NA
## SWL:DRecipientPre NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 15.52 on 54 degrees of freedom
## Multiple R-squared: 0.3019, Adjusted R-squared: 0.2243
## F-statistic: 3.892 on 6 and 54 DF, p-value: 0.002662lm(y ~ D * S)##
## Call:
## lm(formula = y ~ D * S)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.40 12.53 -18.57
## SHH SPS SWL
## 27.17 15.17 NA
## DRecipientPost:SHH DRecipientPre:SHH DRecipientPost:SPS
## -26.56 NA -18.96
## DRecipientPre:SPS DRecipientPost:SWL DRecipientPre:SWL
## NA NA NAlm(y ~ D + S)##
## Call:
## lm(formula = y ~ D + S)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre SHH
## 125.400 9.018 -10.375 5.864
## SPS SWL
## 1.244 NA`?`(glm)
glm(y ~ D)##
## Call: glm(formula = y ~ D)
##
## Coefficients:
## (Intercept) DRecipientPost DRecipientPre
## 125.400 11.441 -8.983
##
## Degrees of Freedom: 60 Total (i.e. Null); 58 Residual
## Null Deviance: 18630
## Residual Deviance: 14470 AIC: 514.7# Description+Sitestart ‘2017_16_07_SSU_Boxplots.R’ Plotting
# richness
plot.df <- ssuguildREADRICH
str(plot.df)## 'data.frame': 183 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
richbox
richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox
richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox
## reads
plot.read <- ssuguildREADRICH
str(plot.read)
reachbox <- ggplot(plot.read, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
reachbox
reachbox <- reachbox + theme_bw(base_size = 15) + xlab("Functional Group") +
ylab("AMF Taxa Reads")
reachbox
reachbox <- reachbox + scale_fill_manual(values = c("red", "darkslateblue",
"gold3"))
reachboxADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA
library(car)
require(MASS)
library(lme4)Choose distribution based on data
find the right distribution for data
ssuguildREADRICHlong$OTU_Richness_Sample.t <- ssuguildREADRICHlong$OTU_Richness_Sample +
1
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "norm")
# LOGNORMAL
qqp(ssuguildREADRICHlong$OTU_Richness_Sample, "lnorm")
# NEG BINOMIAL
nbinom <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Negative Binomial")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "nbinom", size = nbinom$estimate[[1]],
mu = nbinom$estimate[[2]])
# POISSON
poisson <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "Poisson")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "pois", lambda = poisson$estimate)
# GAMMA
gamma <- fitdistr(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma")
qqp(ssuguildREADRICHlong$OTU_Richness_Sample.t, "gamma", shape = gamma$estimate[[1]],
rate = gamma$estimate[[2]])
# modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))
modelingAMF##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH SitePS SiteWL TreatD
## 174.4580 -3.4348 -10.8522 -8.4202 -12.8889
## TreatF TreatO TreatS TreatT
## 11.1111 -24.0580 -0.8703 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 52 Residual
## Null Deviance: 26340
## Residual Deviance: 16940 AIC: 536.3modelingAMF_MP <- glm(OTU_Richness_Sample.t ~ Site + Treat + Description + Rep +
Year, data = ssuguildREADRICHlong, family = gaussian(link = "identity"))
plot(modelingAMF_MP)
stepAIC(modelingAMF_MP)## Start: AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep + Year
##
##
## Step: AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Description + Rep
##
##
## Step: AIC=519.78
## OTU_Richness_Sample.t ~ Site + Treat + Rep
##
## Df Deviance AIC
## - Site 3 12025 519.43
## <none> 10961 519.78
## - Rep 5 16942 536.34
## - Treat 5 20352 547.52
##
## Step: AIC=519.43
## OTU_Richness_Sample.t ~ Treat + Rep
##
## Df Deviance AIC
## <none> 12025 519.43
## - Rep 5 17693 532.98
## - Treat 5 21477 544.81##
## Call: glm(formula = OTU_Richness_Sample.t ~ Treat + Rep, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) TreatD TreatF TreatO TreatS
## 165.889 -12.889 11.111 -28.319 -1.981
## TreatT Rep2 Rep3 Rep4 Rep5
## 1.222 11.737 -8.737 21.431 -1.569
## Rep6
## 35.431
##
## Degrees of Freedom: 60 Total (i.e. Null); 50 Residual
## Null Deviance: 26340
## Residual Deviance: 12020 AIC: 519.4modelingAMF_MP##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat + Description +
## Rep + Year, family = gaussian(link = "identity"), data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH
## 177.106 -6.603
## SitePS SiteWL
## -13.930 -13.110
## TreatD TreatF
## -12.889 11.111
## TreatO TreatS
## -30.959 -2.219
## TreatT DescriptionRecipientPost
## 1.222 NA
## DescriptionRecipientPre Rep2
## NA 11.737
## Rep3 Rep4
## -8.746 21.866
## Rep5 Rep6
## -3.592 39.963
## Year2017
## NA
##
## Degrees of Freedom: 60 Total (i.e. Null); 47 Residual
## Null Deviance: 26340
## Residual Deviance: 10960 AIC: 519.8stepAIC(modelingAMF)## Start: AIC=536.34
## OTU_Richness_Sample.t ~ Site + Treat
##
## Df Deviance AIC
## - Site 3 17693 532.98
## <none> 16942 536.34
## - Treat 5 24729 549.41
##
## Step: AIC=532.98
## OTU_Richness_Sample.t ~ Treat
##
## Df Deviance AIC
## <none> 17693 532.98
## - Treat 5 26339 547.25##
## Call: glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) TreatD TreatF TreatO TreatS
## 166.8889 -12.8889 11.1111 -22.4183 -0.7639
## TreatT
## 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 55 Residual
## Null Deviance: 26340
## Residual Deviance: 17690 AIC: 533modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
data = ssuguildREADRICHlong)
summary(modelingAMFfinal)##
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.980 -8.400 3.265 9.740 32.962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 174.4580 11.3814 15.328 < 2e-16 ***
## SiteHH -3.4348 10.5906 -0.324 0.74699
## SitePS -10.8522 10.2271 -1.061 0.29354
## SiteWL -8.4202 9.9487 -0.846 0.40123
## TreatD -12.8889 8.5090 -1.515 0.13590
## TreatF 11.1111 8.5090 1.306 0.19737
## TreatO -24.0580 8.0233 -2.999 0.00415 **
## TreatS -0.8703 8.7812 -0.099 0.92143
## TreatT 1.2222 8.5090 0.144 0.88634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 325.8149)
##
## Null deviance: 26339 on 60 degrees of freedom
## Residual deviance: 16942 on 52 degrees of freedom
## AIC: 536.34
##
## Number of Fisher Scoring iterations: 2plot(modelingAMFfinal)
modelingAMF <- glm(OTU_Richness_Sample.t ~ Site + Treat, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))
modelingAMF##
## Call: glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) SiteHH SitePS SiteWL TreatD
## 174.4580 -3.4348 -10.8522 -8.4202 -12.8889
## TreatF TreatO TreatS TreatT
## 11.1111 -24.0580 -0.8703 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 52 Residual
## Null Deviance: 26340
## Residual Deviance: 16940 AIC: 536.3stepAIC(modelingAMF)## Start: AIC=536.34
## OTU_Richness_Sample.t ~ Site + Treat
##
## Df Deviance AIC
## - Site 3 17693 532.98
## <none> 16942 536.34
## - Treat 5 24729 549.41
##
## Step: AIC=532.98
## OTU_Richness_Sample.t ~ Treat
##
## Df Deviance AIC
## <none> 17693 532.98
## - Treat 5 26339 547.25##
## Call: glm(formula = OTU_Richness_Sample.t ~ Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Coefficients:
## (Intercept) TreatD TreatF TreatO TreatS
## 166.8889 -12.8889 11.1111 -22.4183 -0.7639
## TreatT
## 1.2222
##
## Degrees of Freedom: 60 Total (i.e. Null); 55 Residual
## Null Deviance: 26340
## Residual Deviance: 17690 AIC: 533modelingAMFfinal <- glm(OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
data = ssuguildREADRICHlong)
summary(modelingAMFfinal)##
## Call:
## glm(formula = OTU_Richness_Sample.t ~ Site + Treat, family = gaussian(link = "identity"),
## data = ssuguildREADRICHlong)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -41.980 -8.400 3.265 9.740 32.962
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 174.4580 11.3814 15.328 < 2e-16 ***
## SiteHH -3.4348 10.5906 -0.324 0.74699
## SitePS -10.8522 10.2271 -1.061 0.29354
## SiteWL -8.4202 9.9487 -0.846 0.40123
## TreatD -12.8889 8.5090 -1.515 0.13590
## TreatF 11.1111 8.5090 1.306 0.19737
## TreatO -24.0580 8.0233 -2.999 0.00415 **
## TreatS -0.8703 8.7812 -0.099 0.92143
## TreatT 1.2222 8.5090 0.144 0.88634
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 325.8149)
##
## Null deviance: 26339 on 60 degrees of freedom
## Residual deviance: 16942 on 52 degrees of freedom
## AIC: 536.34
##
## Number of Fisher Scoring iterations: 2plot(modelingAMFfinal) #Working on examining other variables
modelingAMF <- glm(OTU_Richness_Sample.t ~ Treat + Site, data = ssuguildREADRICHlong,
family = gaussian(link = "identity"))Plotting
plot.df <- ssuguildREADRICH
str(plot.df)## 'data.frame': 183 obs. of 16 variables:
## $ ID : chr "MM933M" "MM933M" "MM933M" "MM936M" ...
## $ Functional.Group : chr "Ancestral" "Edaphophilic" "Rhizophilic" "Ancestral" ...
## $ OTU_Richness_Sample : int 26 2 102 20 1 84 19 1 101 20 ...
## $ Read_Abundance_Sample: num 194 23.2 616.6 103.9 12.2 ...
## $ #SampleID : int 5 5 5 3 3 3 1 1 1 2 ...
## $ BarcodeSequence : chr "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "ATCCCGTATCGATTGG" "GCAACCTTAGAGTGTG" ...
## $ LinkerPrimerSequence : logi NA NA NA NA NA NA ...
## $ Location : chr "E1" "E1" "E1" "C1" ...
## $ Project : chr "MM-Salvage" "MM-Salvage" "MM-Salvage" "MM-Salvage" ...
## $ Year : chr "2015" "2015" "2015" "2015" ...
## $ Site : chr "WL" "WL" "WL" "WL" ...
## $ Treat : Factor w/ 6 levels "C","D","F","O",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ Rep : Factor w/ 6 levels "1","2","3","4",..: 5 5 5 3 3 3 1 1 1 2 ...
## $ HL : num 0.264 0.264 0.264 0.264 0.264 ...
## $ TSDel : chr NA NA NA NA ...
## $ Description : chr "RecipientPre" "RecipientPre" "RecipientPre" "RecipientPre" ...richbox <- ggplot(plot.df, aes(x = Functional.Group, y = OTU_Richness_Sample)) +
geom_boxplot(aes(fill = Description))
richbox
richbox <- richbox + theme_bw(base_size = 15) + xlab("Functional Group") + ylab("AMF Taxa Richness")
richbox
richbox <- richbox + scale_fill_manual(values = c("red", "darkslateblue", "gold3"))
richbox
ADAPTED FROM #GLMM for roots and soil OTU Richness SSU by functional group-modified from “GLM_NAU_6_2017” #—different models for each functional group #Author: Michala Phillips #MODIFY FOR THIS DATA #Tried with negative binomial but it didn’t work
# salvage.div<-adonis2(bird.dist~DIVERSITY, data=birds, permutations = 999,
# method='bray', strata='PLOT') bird.div